% This function runs the second step of the SR approach with specific
% options for what is a allowed to switch per equation
function [theta2,AVarTheta2Step2] = step2SR_threshold_macro_spec(...
    AVarTheta1Step1,outputStep1,bandwidthT,econAct,...
    switch_h0,switch_hx)

% Getting starting values
if any(switch_h0 == 1)
    h0RegimeOn = 1;
else
    h0RegimeOn = 0;
end
if any(any(switch_hx == 1))
    hxRegimeOn = 1;
else
    hxRegimeOn = 0;
end
allowDeSaling  = 1;           % allow for descaling of the var_u to ensure omegeHat is positive definite
thresholdLevel = 0;
[theta2] = Var1Threshold_GMM_closedForm_macro(...
    AVarTheta1Step1.VarFactorsRobust,...
    AVarTheta1Step1.autoCov10Factors,...
    AVarTheta1Step1.autoCov01Factors,...
    outputStep1.factors,outputStep1.macros',...
    econAct,thresholdLevel,allowDeSaling,h0RegimeOn,hxRegimeOn);
% Removing elements in theta2Start to match switch_h0 and switch_hx
nx = length(switch_h0);
for i=1:nx
    if switch_h0(i,1) == 0
        theta2 = rmfield(theta2,['h0',num2str(i),'_2']);
    end
end
for i=1:nx
    for j=1:nx
        if switch_hx(i,j) == 0
            theta2 = rmfield(theta2,['hx',num2str(i),num2str(j),'_2']);
        end
    end
end
selectTheta2 = fieldnames(theta2);
theta2Values = struct2array(theta2)';

% Optimization
if all(switch_h0 ==1) && all(all(switch_hx))
    theta2ValuesOpt = theta2Values;
else
    options         = optimset('Display','iter','MaxIter',5D4,'MaxFunEvals',5D4,'TolFun',1D-10,'TolX',1D-8);
    theta2ValuesOpt = lsqnonlin(@objectFuncStep2Threshold_macro_spec,theta2Values,...
        [],[],options,selectTheta2,AVarTheta1Step1.VarFactorsRobust,...
        AVarTheta1Step1.autoCov10Factors,AVarTheta1Step1.autoCov01Factors,...
        outputStep1.factors,outputStep1.macros',econAct,thresholdLevel,switch_h0,switch_hx);
end
theta2 = values2struct(theta2ValuesOpt,selectTheta2);
% Adding common elements across regimes
for i=1:nx
    if switch_h0(i,1) == 0
        name1 = ['h0',num2str(i),'_1'];
        name2 = ['h0',num2str(i),'_2'];
        theta2.(name2) = theta2.(name1);
    end
end
for i=1:nx
    for j=1:nx
        if switch_hx(i,j) == 0
            name1 = ['hx',num2str(i),num2str(j),'_1'];
            name2 = ['hx',num2str(i),num2str(j),'_2'];
            theta2.(name2) = theta2.(name1);
        end
    end
end

% Ensure that the threshold model is stationary
[delta,theta2] = runInduceStatVAR_threshold_macro_spec(theta2,...
    AVarTheta1Step1.VarFactorsRobust,...
    AVarTheta1Step1.autoCov10Factors,...
    AVarTheta1Step1.autoCov01Factors,...
    outputStep1.factors,outputStep1.macros',econAct,thresholdLevel,switch_h0,switch_hx);

% The standard errors for theta2 using asymptotical inference
AVarTheta2Step2 = getAVarTheta2_threshold_macro_spec(theta2,AVarTheta1Step1.VarFactorsRobust,...
    AVarTheta1Step1.autoCov10Factors,AVarTheta1Step1.autoCov01Factors,...
    outputStep1.factors,outputStep1.macros',econAct,bandwidthT,thresholdLevel,switch_h0,switch_hx);
AVarTheta2Step2.namesFull    = fieldnames(theta2);

% Additional output
AVarTheta2Step2.delta = delta;
AVarTheta2Step2.thresholdLevel = thresholdLevel;
AVarTheta2Step2.switch_h0      = switch_h0;
AVarTheta2Step2.switch_hx      = switch_hx;
end

